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Abstract 

We describe the phase diagram of amorphous solid water by perform- 
ing molecular dynamics simulations using the simple point charge (SPC/E) 
model. Our simulations follow different paths in the phase diagram: isother- 
mal compression/decompression, isochoric cooling/heating and isobaric cool- 
ing/heating. We are able to identify low-density amorphous (LDA), high- 
density amorphous (HDA), and very-high density amorphous (VHDA) ices. 
The density p of these glasses at different pressure P and temperature T agree 
well with experimental values. We also study the radial distribution functions 
of glassy water. In agreement with experiments, we find that LDA, HDA and 
VHDA are characterized by a tetrahedral hydrogen-bonded network and that, 
as compared to LDA, HDA has an extra interstitial molecule between the first 
and second shell. VHDA appears to have two such extra molecules. We ob- 
tain VHDA by isobaric heating of HDA, as in experiment. We also find that 



"other forms" of glassy water can be obtained upon isobaric heating of LDA, 
as well as amorphous ices formed during the transformation of LDA to HDA. 
We argue that these other forms of amorphous ices, as well as VHDA, are not 
altogether new glasses but rather are the result of aging induced by heating. 
Samples of HDA and VHDA with different densities are recovered at nor- 
mal P, showing that there is a continuum of glasses. Furthermore, the two 
ranges of densities of recovered HDA and recovered VHDA overlap at ambi- 
ent P. Our simulations are consistent with the possibility of HDA^LDA and 
VHDA^LDA transformations, reproducing the experimental findings. We 
do not observe a VHDA^HDA transformation, and our final phase diagram 
of glassy water together with equilibrium liquid data suggests that for the 
SPC/E model the VHDA^HDA transformation cannot be observed with the 
present heating rates accessible in simulations. Finally, we discuss the conse- 
quences of our findings for the understanding of the transformation between 
the different amorphous ices and the two hypothesized phases of liquid water. 



I. INTRODUCTION 

Liquids can transform into glass if they are cooled fast enough to avoid crystallization. 
Water is not an exception, although vitrifying liquid water requires very high quenching rates 
(?s 10 6 K/s) [1] — the resulting glass is called hyperquenched glassy water (HGW). Actually, 
water can exist in more than one amorphous form, a property called "polyamorphism" [2-4] . 
Uniaxial compression of ice-I h [5] or ice-/ c [6,7] at temperature T = 77 K to pressures 
P ^ 1 GPa produces a disordered high density material named high-density amorphous ice 
(HDA). If HDA is recovered at T = 77 K and P = 1 bar and then heated isobarically, it 
transforms irreversibly at T ps 125 K to a disordered low-density material named low-density 
amorphous ice (LDA) [5,8]. Both HGW and LDA transform to HDA when compressed at 
T = 77 K at P = 0.6 GPa [8,9]. X-ray and neutron diffraction measurements suggest that 
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LDA is structurally identical to HGW [10] and, although small differences have been found 
[11-13], the common view is that HGW and LDA, are the same material [2,4,14]. 

The LDA-HDA transformation is reversible if it is performed above T m 130 K [15]. 
Furthermore, because the LDA-HDA transformation is very sharp and shows hysteresis, it 
has been suggested that it is a first-order transition [15,16]. It has also been proposed [17] 
to interpret the LDA-HDA transition as the glass counterpart of the first order transition 
between two liquids of different structure, the liquid-liquid phase transition scenario [18]. 
In this view, LDA is the glass associated with a low-density liquid (LDL) while HDA is the 
glass associated with a high-density liquid (HDL) [16,17,19-21]. 

Computer simulations have been able to reproduce many of the experimental results. 
More than 15 years ago, simulations using the TIP4P model found a transformation of ice- 
i/j into HDA when compressed at T = 80 K and a transformation of HDA to LDA upon 
isobaric heating at zero P [22]. The LDA-HDA reversible transformation is also found in 
simulations using the ST2 [18] and TIP4P [17] models. Computer simulations in the liquid 
phase also show a LDL-HDL first-order transition [18,23-25]. 

Very recently, a glass denser than HDA was identified in experiments, and named very- 
high density amorphous (VHDA) ice [26]. VHDA is obtained by isobaric heating at high 
P of HDA from T = 77 K to T 165 K. Two samples of HDA, one heated isobarically at 
P = 1.1 GPa and the other at P = 1.9 GPa, both relax at T = 77 K and P = 1 bar to 
the same structure with a density of p = 1.25 ± 0.01 g/cm 3 , i.e., 9% denser than HDA 
(p = 1.15 g/cm 3 ) and pa 40% denser than LDA (p = 0.94 g/cm 3 ). It appears that the 
properties of HDL are closer to VHDA than to HDA [27,28]. In recent reports [29,30], it 
has been shown that VHDA can be obtained in computer simulations by following the same 
procedure as in experiments. Further, we argued [30] that VHDA may not be a new glass 
different than HDA, but rather VHDA (and not HDA) results from partial annealing of the 
"more metastable" HDA structure. Further, our simulations also showed that VHDA is the 
glass obtained by isobaric quenching of the equilibrium liquid at high P. 

It is an open question whether all forms of glassy water can be classified into two amor- 
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phous forms, which we call LDA and VHDA. While this might be the case for glasses 
obtained by cooling liquid water or compressing ice, there are many different experimental 
techniques to obtain glassy water. For instance, amorphous ice can be also obtained by 
exposing crystalline ice to radiations such as electrons [31], ultraviolet photons [32], and ion 
bombardment [33]. The situation is more complicated when considering the effect of aging 
in the glassy state. Recent works [34-36] suggest that aging takes place on the experimental 
time scale when HDA is kept at ambient P and T is below the HDA-LDA transformation 
T. The p and structure of the sample changes with the annealing T and annealing time. 

In this work we rebuild the phase diagram of amorphous water on the basis of simula- 
tions using the SPC/E model. Our simulations consist mainly of isobaric cooling/heating 
and isothermal compression/decompression cycles. In section II, we describe briefly the 
simulation details. In section III, we show that compression of LDA leads to HDA. Isobaric 
heating at different Ps is reported in section IV, where the HDA^VHDA transformation is 
also described. In section V, we describe the results of recovering HDA and VHDA at normal 
P. The structures of all glasses are studied in section VI. The HDA^LDA, VHDA^LDA 
and VHDA— >HDA transitions are investigated in section VII. Finally, in section VIII we 
discuss our results and recent experiments in the context of a possible first order transition 
line separating two amorphous ices at low T, and two liquid phases at high T. 

II. SIMULATIONS 

We perform molecular dynamics (MD) simulations using the extended simple point 
charge (SPC/E) model of water [37]. This model has been extensively used to study the 
thermodynamics [38,39] and dynamics [40,41] of liquid water. Thermodynamic and dynam- 
ical properties of SPC/E water are well known, and are consistent with experimental facts. 
The SPC/E model reproduces the thermodynamic anomalies characterizing water (e.g. it 
shows a maximum in p [42]). This model has also been used to study glassy water, and 
is able to produce glassy states corresponding to LDA, HDA and VHDA [30,43]. Here, we 
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simulate a system oi N = 216 molecules [44], and average the results over 16 independent 
trajectories. We use periodic boundary conditions. Long range forces are treated using the 
reaction field method. 

To study the phase diagram of glassy water we perform three different kinds of sim- 
ulations: (i) isobaric cooling/heating, (ii) isochoric cooling/heating, and (iii) isothermal 
compression/decompression. 

During cooling/heating at both constant p or P, we change T by rescaling the velocities 
of the molecules. At every time step St — 1 fs, we increase the T by 5T = q5t, where q 
is the cooling/heating rate and q = ±3 x 10 10 K/s. These values for q have been used in 
previous simulations to study the effect of cooling and heating rates on the glass transition 
temperature T g = T g (P) [45,46]. The glass obtained with a cooling rate q = — 3 x 10 10 K/s 
behaves, upon heating at q = +3 x 10 10 K/s, as a slow-cooled glass as found in experiments 
showing no signs of hyper quenching effects. On the other hand, the lowest cooling/heating 
rates accessible nowadays in computer simulations are \q\ as 10 10 K/s [28]. 

For the compression/decompression simulations at T > K we perform MD simulations 
at constant p for intervals of 1 ps. At the end of each interval of 1 ps, we increase p by 5p = 
±5 x 10" 5 g/cm 3 ,so our compression/decompression rate is 5p/5t = ±5 x 1CT 5 g/cm 3 /ps. 
This value was already used to study the potential energy landscape for the LDA-HDA 
transformation [43]. These p changes are performed by rescaling isotropically the coordinates 
of the molecular center of mass. For the compression simulations at T = K, at each 
simulation step we change p by 5p = ±5 x 10~ 5 g/cm 3 and then minimize the energy. At 
each step p is modified by rescaling isotropically the center of mass of each molecule. 

III. LDA-HDA TRANSFORMATION 

We prepare LDA by hyperquenching at constant p = 0.9 g/cm 3 equilibrium liquid con- 
figurations obtained at T = 220 K. We have chosen a glass obtained quenching a low T 
liquid equilibrium configuration — i.e. properly speaking a HGW — -at p = 0.9 g/cm 3 as 
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starting LDA structure to guarantee that the chosen configuration is a low p glass with an 
optimized tetrahedral network of hydrogen bonds. To obtain HDA, we compress LDA at 
constant T [43] from p = 0.9 g/cm 3 to p > 1.5 g/cm 3 . The T is fixed at T = 0, 77, 100, 
120, 140, 160, and 170 K. These T are below T g (P). Evaluation of the specific heat upon 
heating the glass with a heating rate of q = +3 x 10 10 K/s to the liquid phase indicates that 
T g (P) > 170 K for all Ps studied. 

In Fig. 1 we show P(p) during isothermal compression at four of the simulated T. We 
also show equilibrium data taken from the liquid phase simulations reported in Ref. [39]. 
From Fig. 1 we note: 

(i) For all T < 170K isotherms show a transformation from LDA to HDA. The linear 
increase of P(p) crosses to a transition region around p c ps 1.03 g/cm 3 and P c ps 
0.7 GPa for T = K and p c ps 0.97 g/cm 3 and P c ps 0.13 GPa for T = 170K. The 
LDA-HDA transition starts at higher p and P on decreasing T. This feature has been 
observed in computer simulations of different models of water, including ST2 [18,17] 
and TIP4P [17], and agrees with experiments [15,47]. 

(ii) All the compression curves at T > K collapse to a single curve at very high p, 
suggesting that there is a single HDA state at very high P. This is not the case 
for LDA because the P of LDA while compressing at p = 0.9 g/cm 3 depends on 
T. Furthermore, a study of the potential energy landscape during the LDA-HDA 
transformation shows clear differences between two LDA samples at low P and T [43]. 

(iii) The isotherms in the glassy states cross each other at approximately p = 0.98 g/cm 3 
and P = 0.2 GPa (see also Fig. 7). This could suggest (if these out of equilibrium 
isotherms could be interpreted as an effective equation of state of the glass) that p 
increases upon isobaric heating the glass at high P, while p decreases upon isobaric 
heating the glass at low P. 

(iv) We note that as the compression temperature increases from T = 0KtoT = 170 K the 
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isotherms in the glassy state approach the liquid isotherm for T = 210 K. Accordingly, 
we see that at low p (LDA, e.g. p = 0.9 g/cm 3 ) increasing T increases P, while at 
high p (HDA, e.g. at p = 1.3 g/cm 3 ) increasing T decreases P. 

IV. ISOBARIC HEATING OF GLASSY WATER. VERY HIGH DENSITY 
AMORPHOUS ICE AND "OTHER" AMORPHOUS ICES. 

A. Heating at High P 

Experimentally VHDA is obtained by isobaric heating of HDA (T = 77 K, p = 
1.33 g/cm 3 ) at P = 1.1 GPa and HDA (T = 77 K, p = 1.45 g/cm 3 ) at P = 1.9 GPa 
[26]. The products are VHDA (T w 160 K, p w 1.37 g/cm 3 ) and VHDA (T w 177 K, 
p ps 1.51 g/cm 3 ), respectively. Simulations with the SPC/E model also show a transforma- 
tion between HDA to another glassy state identified with VHDA [30]. These simulations 
also suggest that HDA and VHDA are not different glasses, but rather that VHDA is the 
result of the relaxation of HDA upon heating. Furthermore, it is found that VHDA is the 
same glass obtained by cooling HDL [28,30]. 

We obtain HDA by compressing LDA at T = 77 K as explained in the last section. We 
confirm the transformation of HDA to a denser glass by isobaric heating of HDA at different 
Ps. The HDA configurations correspond to four different states, indicated by state points 
H, I, J and K in Fig. 2 (see Table I). For the four cases studied the densities upon heating 
are shown in Fig. 3. In all cases, p increases when heating the glass from T = 77 K to 
T ~ 170 K still below T g (P). On increasing T, the glass melts to the equilibrium liquid. 
The p increase is smaller the higher the P and it is indicated by arrows in Fig. 2. Note that 
the increase in p of ~ 0.04 g/cm 3 that we find at P = 1.38 GPa is very close to the estimated 
experimental values [26] of 0.04 g/cm 3 and 0.06 g/cm 3 for isobaric heating at P = 1.1 and 
P = 1.9 GPa, respectively. 

The p(T) values shown in Fig. 3 at high T are the ones obtained in the equilibrium liquid. 
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Indeed, at P = 1.38 GPa and T = 230 K the relaxation time r is r « 100 ps [48]. Heating 
at a rate of g fc = 3 x 10 10 K/s, the T-change during 100 ps is about 3 K. Therefore, the 
system has time to relax before T changes considerably and the p measured is the same as 
in the equilibrium liquid. The decrease of p observed in Fig. 3 for T > 230 K (in the liquid 
phase) is in agreement with the location of the isotherms for T = 210 K and T = 230 K at 
high P shown in Fig. 2. 

Computer simulations show that the final glass connected to the liquid phase by cooling 
at high P is VHDA [30]. We argue [30] that HDA is less stable than VHDA and it is a 
kinetically trapped glass that could not relax during the compression of LDA. The hypothesis 
that HDA is a kinetically trapped glass was already proposed based on recent computer 
simulations performed with the TIP4P model [29] and it was first proposed by Mishima 
who observed that HDA samples (made by different P-T conditions) heated or annealed 
to 130 — 150 K at 1 — 1.5 GPa were characterized by identical X-ray patterns [49]. To 
provide further support to this hypothesis, we age configurations corresponding to state 
point H in Fig. 2 before heating, i.e. at P = 1.38 GPa and T = 77 K. The evolution of 
p with time t is shown in Fig. 4. As t goes on, p increases approaching the value of p 
corresponding to VHDA at P = 1.38 GPa. Another support to the hypothesis that VHDA 
can be interpreted as 'relaxed HDA' comes from the study of the effect of the compression T 
in the LDA— >HDA transformation. If HDA is less stable than VHDA, we should find that 
(for a given final P after compression of LDA) the larger the compression T, the closer is 
HDA to the corresponding VHDA. In other words, the higher the compression T of LDA, 
the more the system can relax to a more stable state. Accordingly, Fig. 2 shows that as the 
compression T is increased, the compression isotherms starting from LDA shift (for a given 
P) to higher p, approaching VHDA. 
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B. Heating at Intermediate P 



We explore now the effect of isobaric heating of the glass at intermediate P, in the LDA- 
HDA transition region, indicated by state points D, E, F, and G in Fig. 5 (see Table I). The 
lower the P of these HDA forms, the more contaminated they are by LDA. Upon isobaric 
heating from T = 77 K up to T ps 170 K, the states indicated by state points D, E, F, 
and G shift to higher densities as indicated by the arrows in Fig. 5. The evolution of the 
p upon heating these four glasses are shown in Fig. 6. The situation is analogous to the 
HDA-VHDA transformation discussed above. The increase in T allows these glasses to relax 
to more stable glass configurations characterized by a higher p. The final glasses found after 
heating at these intermediate P correspond to the "VHDA" obtained at high Ps. We note 
also that at a given intermediate P (e.g. P = 0.8 GPa), the isotherms obtained in the LDA- 
HDA transformations shift to higher p as the compression T increases. This is because, 
as discussed in Sec. IV A, the higher the compression T, the more the glass relaxes during 
compression. 



C. Heating at Low and Negative P 

Isobaric heating at Ps below the P at which the glass isotherms cross (P ~ 0.2 GPa, 
p ps 0.98 g/cm 3 ) presents some differences from the cases studied above. We discuss the 
heating from three LDA configurations indicated by state points A, B and C in Fig. 7 (see 
Table I). The heating of LDA represented by state point C at P = 0.01 GPa is shown 
in Fig. 8(a). We observe a decrease in p from 0.96 g/cm 3 at 77 K to p ps 0.952 g/cm 3 
at T = 150 K, and then an increase in p until T = 230 K (already in the liquid phase). 
This behavior of LDA upon heating at P = 0.01 GPa is consistent with the location of 
the isotherms in the glassy state. Isotherms in the glassy state shift to lower densities 
as T increases, suggesting that the more stable glasses at P = 0.01 GPa are at lower 
densities. Therefore, the glass represented by state point C moves upon heating first to 
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lower p approaching the glass isotherm at T = 170 K in the glassy state, and then shifts 
back to higher p approaching the T = 230 K equilibrium liquid isotherm. The relaxed 
glass obtained before the glass transition is a "very-low density glass" ("VLDA") and is 
conceptually analogous to the VHDA obtained at high P. A more precise nomenclature for 
these glasses might be relaxed-LDA (RLDA) and relaxed-HDA (RHDA), respectively (see 
also [35]). 

A situation similar to the one discussed previously for state point C in Fig. 7 holds 
for state point B in Fig. 7, i.e. for the isobaric heating at P = —0.17 GPa of LDA (p = 
0.94 g/cm 3 ). The p upon heating for this case is shown in Fig. 8(b). We observe an initial 
decrease in p down to xs 0.93 g/cm 3 at T « 170 K approaching the value of p ps 0.92 g/cm 3 
corresponding to the glass isotherm at T = 170 K. As T increases up to 230 K, a non 
monotonic density dependence arise, probably because the p of the equilibrium liquid at 
P = -0.17 GPa and 230 K is also m 0.93 g/cm 3 . The final decrease in p for higher T is 
due to the fact that the liquid isotherms at P = —0.17 GPa in Fig. 7 shift to lower p as T 
increases. 

Finally, Fig. 8(c) shows the effect of heating LDA at P = —0.55 GPa, indicated by state 
point A in Fig. 7. In this case, p decreases abruptly upon heating and LDA transforms into 
gas (i.e., it "sublimates") because there are no liquid isotherms at this P. 

V. RECOVERING OF HDA AND VHDA AT AMBIENT P. A CONTINUUM OF 

GLASSES. 

In this section, we investigate the effects of isothermal expansion of HDA and VHDA 
at low T down to p = 0.85 g/cm 3 . The structure of HDA and VHDA has been described 
at T = 77 K and P = 0.01 GPa [26,50]. To compare our results with the experiments, we 
also recover HDA and VHDA at these conditions. For both HDA and VHDA, the system is 
decompressed at constant T = 77 K (in the case of VHDA, we first bring the system back 
to T = 77 K by keeping the P constant). We recover the configurations of HDA represented 
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by state points H, I, J, and K in Fig. 2, and the corresponding configurations of VHDA. 

We first discuss state points H and I, the cases of HDA at P = 1.38 GPa and P = 1.9 GPa 
[see inset of Fig. 9(a)]. State points H' and F are the glasses obtained after heating HDA, 
i.e., VHDA 165 K). In agreement with experiments, cooling VHDA back to T = 77 K 
increases p. The final glasses after cooling, i.e., VHDA 77 K), are indicated by state 
points H" and I" in the figure. We note that VHDA(^ 77 K) can be cycled back upon 
heating to VHDA 165 K). This is also in agreement with the experimental results [26]. 
Figure 9(a) shows the evolution of P as a function of p during decompression of the system 
starting at state points H, I (for HDA), H" and I" (for VHDA). During decompression, 
VHDA is denser than the corresponding HDA, and only at P xs —0.15 GPa the densities of 
the decompressed samples of VHDA are equal. 

In experiments, HDA and VHDA recovered at ambient P and T = 77 K are characterized 
by only one state with p = 1.17 g/cm 3 and p = 1.25 g/cm 3 , respectively. In Fig. 9(b) 
and 9(c) we show the P during decompression of the four HDA and the four VHDA forms, 
respectively. The recovered HDA at P = GPa and T = 77 K have different densities, 
and do not collapse to a single state. Also the recovered VHDA forms are characterized by 
different densities, notwithstanding the same P and T. We observe that at P GPa, 
HDA densities fall in the interval 1.15 — 1.24 g/cm 3 while VHDA densities fall in the interval 
1.22 — 1.28 g/cm 3 . Furthermore, these two p intervals overlap, meaning that recovered HDA 
and VHDA are indistinguishable by their p at zero P. It is also clear from our results that 
there is a continuum of glasses at normal P. A continuum of glasses at normal P has been 
already suggested by experiments in Refs. [34,35] and obtained in computer simulations 
using the TIP4P model [29]. Our simulations also show that as P decreases to negative 
values, all the recovered HDA and VHDA forms collapse to a single state at p = 1.05 g/cm 3 
and P ~ -0.4 GPa. 
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VI. RADIAL DISTRIBUTION FUNCTIONS 
A. LDA and HDA at normal P 

In this section, we study in detail the radial distribution function (RDF) of amorphous 
ice. We compare the results from simulations and experiments on LDA, HDA and VHDA. 
We also investigate how the structure of HDA changes as HDA is transformed to VHDA at 
high P. The RDFs for glassy water are commonly measured at T xs 80 K and ambient P. 
In Fig. 10(a) we show the RDF obtained from simulations for LDA (p = 0.96 g/cm 3 ) and 
recovered HDA (p = 1.2 g/cm 3 ) both at P = 0.02 ± 0.02 GPa and T = 77 K. Recovered 
HDA is obtained by isothermal decompression at T = 77 K of HDA (p = 1.37 g/cm 3 ) [state 
point H in Fig. 9(a)], as discussed in the section above. Figure 10(a) can be compared 
with Fig. 1 in Ref. [50], which shows the experimental RDF of LDA and recovered HDA at 
the state point (P = 1 atm, T = 80 K). Both figures are very similar. For example, the 
maxima of goo{ r ) for LDA in our simulations are located at 0.27, 0.44 and 0.66 nm, while 
experimentally [50] they are located at 0.27, 0.45 and 0.68 nm. When going from LDA to 
HDA by compression, the first peak of goo( r ) increases, the second peak gets wider and 
smaller, and the third peak shifts to lower values of r. However, for HDA we find some 
differences between experiments and our simulations. The second peak of goo( r ) for our 
HDA simulations cannot be resolved into two peaks as clearly as in experiments [50], and 
the first minimum in goo(r) is located at r = 0.30 nm instead of at ~ 0.32 nm. The same 
observations with respect to the second peak in goo( r ) can be made from MD simulation 
using the TIP4P model [22,29] and from a reverse Monte Carlo simulation [51]. We note that 
the second peak of goo{ r ) is very sensitive to T; measurements at P = GPa, T = 100 K, 
p = 1.16 g/cm 3 show that the second peak of goo( r ) cannot be resolved into two peaks [52]. 

The plots for gou(r) for LDA and HDA in Figure 10(a) show the same features as those 
found experimentally in Ref. [50]. When going from LDA to HDA, the first peak of gon(r) 
increases, the second peak gets wider, and the third peak practically disappears. However, 



12 



we note that for both LDA and HDA the relative heights of the first and second peaks 
of gou( r ) in experiments and simulations are different. Experiments [50] show that both 
peaks have almost the same amplitude for LDA, while for HDA the first peak is smaller 
than the second one. Instead, SPC/E shows a larger first peak in <?oh(?") for both LDA and 
HDA suggesting that the structure is more tetrahedral than in real water. In accord with 
experiment [50], for the case of gnn(r) we find an increase of the first peak and decrease of 
the first minimum of the RDF. 

B. HDA and VHDA at ambient P 

Figure 10(b) shows the RDF of recovered HDA (p m 1.37 g/cm 3 ) and recovered VHDA 
(p = 1.28 g/cm 3 ) generated as discussed in the previous sections (P = 0.02±0.02, T = 77 K). 
VHDA was obtained by isobaric heating of HDA (p ps 1.37g/cm 3 ) up to T ps 165 K. 
Figure 10(b) can be compared with Fig. 1 in Ref. [27], which shows the experimental RDF 
for recovered HDA and VHDA at (P = 1 atm, T = 80 K). In accord with experiments, 
when going from HDA to VHDA, goo (r) shows a weak decrease in the first peak and a 
shift of the third peak to lower values of r. Furthermore, the second wide peak of goo( r ) 
for HDA practically disappears at r ps 0.43 nm. In the experimental goo( r ), the second 
peak in HDA merges into the first peak, filling the gap at r = 0.33 nm, the first minimum 
for HDA. As a consequence, the first peak develops a shoulder with a second maximum at 
r = 0.35 nm. In our simulations, this does not occur. Instead, the second maximum is 
located at r « 0.33 — 0.34 nm and the first minimum in HDA is not filled, i.e., we still 
can see a clear minimum in goo(r) for VHDA. The presence of such a minimum in the 
goo(r) for VHDA has been observed also in MD simulations using the TIP4P model [29]. 
When moving from HDA to VHDA, we find that gon(r) barely changes. There is a weak 
shifting of the second peak toward lower values of r as it is found experimentally [27], and 
the observed decrease of the first peak is barely seen in our simulations. We find that gmi(r) 
also shows weak changes. As in [27], we find a decrease in the first peak and increase of the 
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first minimum upon heating. 

We expect to find differences in the RDF between experiments and simulations. Prop- 
erties for the SPC/E potential at a given T correspond to those in experiments found at 
approximately T' T + 40 K [18]. As we already mentioned, the goo(r) we find for HDA 
at T = 77 K is more similar to that found in experiments at T = 100 K [52] than the 
corresponding one at T = 80 K [50]. Recovering VHDA at higher T will allow neighbor 
oxygen atoms to come closer, thereby filling the gap at r ss 0.3 nm in goo( r )- Therefore, 
differences between simulations and experiments in the RDF of VHDA are expected to be 
smaller at higher T. 

C. Isobaric Heating of HDA to VHDA at high-P 

In Fig. 11 we show the evolution of the RDF on heating HDA (p = 1.37 g/cm 3 ) at P = 
1.9 GPa from T = 77 K up to T ps 175 K. Figure 11(a) shows goo( r )- Upon heating HDA, the 
first peak in goo( r ) decreases and molecules move, filling the gap at r = 0.29 nm. The second 
peak shifts to lower values of r, merging with the first peak, and the minimum at r = 0.41 nm 
decreases. This RDF of VHDA is remarkably similar to the experimental distributions [27] 
for recovered VHDA at P = 0.02 ± 0.02 and T = 77 K, suggesting that differences observed 
in our recovered configurations set up during the decompression process. Figure 11(b) 
shows the evolution of gon(r) upon heating. The distribution shows a decreasing first 
peak and increasing first minimum. These changes could not be observed in the recovered 
configurations of HDA and VHDA at P = 0.02 GPa. In Fig. 11(c) we note that gnn{r) 
presents mainly a weak increase of the first minimum. 

D. Coordination Number of LDA, HDA and VHDA 

Integration of the experimental gon(r) at T = 80 K and normal P for LDA, HDA, and 
VHDA in the range 0.14 nm < r < 2.5 nm indicates that each oxygen (O) atom is surrounded 
by 2 ± 0.1 hydrogen atoms [27,50]. Furthermore, integration of goo( r ) for r < 3.1 nm (the 
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location of the first minimum in the experimental VHDA goo( r )) indicates that each O has 
4 nearest-neighbor O atoms. Therefore, LDA, HDA and VHDA are characterized by a 
tetrahedrally-coordinated fully hydrogen bonded network. 

Our simulations agree with these experimental findings: we find from Fig. 10 that at 
P = 0.02 GPa and T = 77 K in LDA, HDA and VHDA each oxygen is surrounded by 
1.95 ± 0.5 hydrogen atoms (for r < 2.5 nm) and by 4 — 4.1 oxygen atoms for r < 0.303 nm 
(the first minimum of the simulated goo(r) for VHDA). Experiments [27,50] also indicate 
that differences in structure of HDA and VHDA as compared to LDA basically arise from 
the presence of extra interstitial neighbors. Integration of the experimental goo{ r ) from 
r = 0.31 nm (the first minimum of goo{ r ) f° r VHDA) to r = 0.33 nm (the first minimum of 
goo( r ) fc> r HDA), indicates that, in this range of distances, an O atom has 0.9 1) oxygen 
neighbors in HDA and 1.7 (~ 2) in VHDA. In simulations, integrating from r > 0.303 nm 
(first minimum of goo(r) for VHDA) up to r = 0.33 nm, we find that each O atom has 1.1 
(« 1) neighbor O in HDA and 1.8 (« 2) in VHDA (in LDA this value is only 0.1). Therefore, 
simulations are consistent with experiments showing that, in comparison to LDA, there is 
one extra interstitial molecule in the structure of HDA and approximatively two in the 
structure of VHDA. 

VII. THE HDA-LDA, VHDA-LDA AND VHDA-HDA TRANSFORMATIONS 

A. HDA-LDA transformation 

It was first observed by Mishima et al. [5] that recovered HDA (p xs 1.17 ± 0.02 g/cm 3 ) 
at T = 77 K and ambient P transforms to LDA (p pa 0.94 ± 0.02g/cm 3 ) upon isobaric 
heating to T ~ 117 K. In Ref. [22] it is found that simulations at ambient P also show 
this transition between HDA and LDA. Here we show that our simulations at ambient P 
show remarkably similar results to those found in [22]. However, we point out that a direct 
HDA-LDA transition with the SPC/E model at the present heating rate qh = 3 x 10 10 K/s 
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takes place only at negative P. 

Configurations of HDA (p = 1.4 g/cm 3 ) are obtained by isothermal compression of LDA 
at 77 K up to P = 2.14 GPa with a rate of dp/dt = 5 x 10~ 5 g/cm 3 /ps, as explained in 
Sec. II. HDA is recovered by isothermal decompression at T = 77 K down to P —0.7 GPa 
with a rate of dp/dt = — 5 x 1CT 5 g/cm 3 /ps [see Fig. 12(a)]. We then heat isobarically 
recovered HDA (p = 1.2 g/cm 3 ) at P = 0.01 GPa. The evolution of the system in the 
P-p plane is indicated by the upper horizontal arrow in Fig. 12(a). p{T) upon heating at 
P = 0.01 GPa is shown in Fig. 12(b). We observe in Fig. 12(b) that p « 1.00 g/cm 3 at 
T = 220 K, while the p of LDA (T = 77 K, P = 0.01 GPa) is w 0.96 g/cm 3 . Therefore, 
with the present rate ^ = 3 x 10 10 K/s, even when the system becomes less dense upon 
heating, it does not reach the p corresponding to LDA at 77 K and 0.01 GPa. Simulations 
using the TIP4P model are similar (Fig. 1 in Ref. [22]). In Ref. [22], the p of the system at 
T 220 K (during the heating process) is even larger than 1 g/cm 3 . The p corresponding 
to LDA in Ref. [22] occurs when the T is T > 260 K and the system is in the liquid phase. 

To find the HDA-LDA transition at P = 0.01 GPa, the system has to "move" in Fig. 12(a) 
(upper horizontal arrow) from p = 1.2 g/cm 3 to p = 0.96 g/cm 3 crossing liquid isotherms, 
like the 210 K-liquid isotherm shown in the figure. In the previous paragraph, we show 
that at the present qh the system gets trapped at a state with p corresponding to the liquid 
phase at T = 210 K. It is not clear whether a much slower q h could allow the system to 
reach the p of LDA. To find a direct HDA^LDA transformation with q h = 3 x 10 10 K/s 
we heat recovered HDA at negative Ps. This is motivated because at P < the system 
should cross first the glass isotherms and then the liquid isotherms to reach the liquid phase. 
We confirm this view by heating isobarically recovered HDA at P = —0.3 GPa [see lower 
horizontal arrow in Fig. 12(a)]. p(T) upon heating is shown in Fig. 12(c) and is consistent 
with the possibility of a HDA— >LDA transformation at P = —0.3 GPa, without any liquid 
intermediate state. 

We note that we also repeated the isobaric heating procedure at P = —0.71 GPa, where 
"recovered HDA" (which is already LDA) has a p of 0.82 g/cm 3 . On isobaric heating p 
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decreases abruptly and the glass transforms to gas. In this case, the system on the P-p 
plane moves away from the LDA isotherms and the glass eventually transforms into gas 
because there are no liquid isotherms at this P. 

B. VHDA-LDA transformation 

Experiments show that isobaric heating of "recovered VHDA" at P = 0.11 GPa from 
T = 77 K to T Rs 127 K produces LDA [26]. We test this experimental finding with our 
simulations by heating configurations of VHDA at P = 0.11 GPa. We generate VHDA 
configurations by heating HDA (p = 1.31 g/cm 3 , P = 1.38 GPa) from 77 K to T « 160 K. 
We then recover VHDA at T = 77 K, P = 0.11 GPa and p = 1.27 g/cm 3 . 

Upon isobaric heating at P = 0.11 GPa (see upper horizontal solid arrow in Fig. 13(a)), 
the system expands as shown in Fig. 13(b). However, the system does not reach the LDA 
state. For this, the final p before the glass transition T should be 0.97 g/cm 3 , the p of 
LDA [see Fig. 13(a)]. Instead, Fig. 13(b) shows that the p at T « 170 K is « 1.12 g/cm 3 , 
close to the p of equilibrium liquid at T = 210 K and P = 0.11 GPa — suggesting that 
the isobaric heating of the recovered VHDA in our simulations brings the system to the 
liquid state without an intermediate LDA. The situation is similar to the above discussion 
of the HDA— >LDA transformation at P = 0.11 GPa. The glass upon heating evolves to the 
liquid phase with no intermediate transformation to LDA. In analogy with the situation of 
the HDA^LDA transition, a VHDA— >LDA transition could be obtained, with the present 
heating rates, at negative Ps, i.e., at P such that the system crosses first the glass isotherms 
and then the liquid isotherms. The p upon isobaric heating at P = —0.3 GPa of recovered 
VHDA is shown in Fig. 13(c) [see also the lower horizontal arrow in Fig. 13(a)] and it is 
consistent with the possibility of a VHDA^LDA transformation. 
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C. VHDA-HDA transformation 



Experiments show that isochoric heating of recovered VHDA at P = 0.02 GPa and 
T = 77 K up to T = 140 K appears to produce HDA [26]. In fact, the x-ray pattern of the 
final sample has been interpreted as a mixture of VHDA and HDA. We use our simulations to 
test for the presence of VHDA— >HDA transformation by heating configurations of VHDA at 
constant volume. Configurations of VHDA are obtained by isobaric heating of HDA (77 K, 
1.31 g/cm 3 ) at P = 1.38 GPa up to T w 160 K. VHDA is then recovered at (P = 0.02 GPa, 
T = 77 K, p = 1.26 g/cm 3 ). The evolution of the system in the P — p plane during isochoric 
heating at p = 1.26 g/cm 3 from T = 77 K up to T ~ 170 K is indicated by the vertical arrow 
in Fig. 13(a). P(T) is shown in Fig. 13(d). P increases monotonically upon heating and 
approaches the value 0.72 GPa at T 200 K. At this P and T, the same thermodynamic 
properties as the equilibrium liquid at T = 210 K are recovered, as indicated by the liquid 
isotherm in Fig. 13(a). Therefore, our simulations do not show a VHDA-HDA transformation 
because the system melts in the liquid phase. The situation is similar to the HDA^LDA 
transformation at P = 0.01 GPa discussed above. The P increase during isochoric heating 
is a consequence of the VHDA^liquid transition. 

It is not clear whether a much slower heating rate (nowadays, not accessible in computer 
simulations [28]) will show a VHDA-HDA transformation or whether such a transformation 
does not exist. In the short term, experiments may help clarify this issue. 

VIII. DISCUSSION 

The present study was inspired by several motivations. Among them were (i) the goal of 
testing how closely simulations based on the SPC/E potential reproduce the experiments, 
and (ii) the goal of deriving a concise view of the phase diagram of amorphous ice. 

We find that even using the simple SPC/E model for water, simulations of isothermal 
compression of LDA produce another glass, HDA. By repeating the process at different 
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compression T we find that the glass isotherms cross one another in the P-p plane. This 
leads to the different behavior of glassy water at low and high P. 

We investigate how under isobaric heating HDA converts into a denser glass, VHDA, 
while under isobaric heating LDA evolves toward a less dense state. The higher is the 
compression T of LDA, the higher is the p of the final HDA form at a given P, and the 
closer it is to the corresponding VHDA. This finding is consistent with the view of VHDA 
as the result of relaxation of HDA upon isobaric heating — i.e., the higher the compression 
T, the more the glass is able to relax during compression, approaching VHDA. In fact, by 
aging HDA at high P before heating to obtain VHDA we observe an increase of p with time 
suggesting that VHDA is the result of relaxation of HDA. A similar interpretation of HDA 
as an unrelaxed glass was proposed by Mishima when forming HDA by compression of ice 
I h at different T [49]. 

We also investigate the thermal stability of the glasses formed during the compression 
process of LDA. To do this, we repeat the protocol to transform HDA into VHDA at high-P. 
In other words, we perform isobaric heatings of samples of LDA and samples obtained during 
the transformation of LDA to HDA. As T increases, irreversible relaxation processes take 
place generating annealed versions of the starting samples. The resulting annealed glasses 
(obtained at low and intermediate P) are analogous to the VHDA (obtained at high P). 

Densities at different Ps and Ts for LDA, HDA and VHDA agree well with experiments. 
Furthermore, comparison of the structure of these glasses at the level of the radial distribu- 
tion functions shows differences between the three forms — LDA, HDA, and VHDA — similar 
to those observed in experiments. In particular, we find that LDA, HDA and VHDA are 
characterized by a tetrahedral hydrogen bonded network. As compared to LDA, HDA has 
an extra interstitial neighbor between the first and second shell while VHDA has two extra 
molecules. 

Stimulated by experimental work, we also simulate the recovering process of different 
samples of HDA and VHDA at normal P. We find a continuum range of densities for 
each family of glasses. Furthermore, the ranges of densities of recovered HDA and recovered 
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VHDA overlap at normal P. A continuum of glasses was proposed in Ref. [34] and mentioned 
as a possibility in Ref. [26] . 

We discuss the LDA, HDA and VHDA results with respect to the equation of state 
of liquid water. The comparison of equilibrium liquid isotherms with those for the glassy 
state in our simulations allow one to predict the effect of cooling/heating and compres- 
sion/decompression of amorphous and/or liquid water. For example, one can understand 
why it is possible to find HDA^LDA and VHDA^LDA transformations in SPC/E at 
P = —0.3 GPa and why it is not the case at P = 0.01 GPa (at least with the present 
heating rate). We note that experimentally these two transformations to LDA occur at 
P ks GPa. However, we expect shifts in the values of P, T and p when comparing experi- 
ments and simulations. The possible VHDA^HDA transformation [26] is not observable in 
simulations with the SPC/E model with the present heating rate. 

Some experiments and computer simulations performed in the last few years have chal- 
lenged the phase diagram proposed in [4,17,18]. Aging effects and the previous history of 
the glass (especially at T < T g {P)) make it difficult to characterize the glass state in an 
unique way. In fact, a glass depends on how it is prepared: cooling/heating rates, com- 
pression/decompression rates and aging times are variables that have to be included when 
describing the glass state [46,53]. Therefore, a phase diagram for a metastable system has to 
be viewed with caution. In the phase diagram proposed in [4,17,18], a first order transition 
line separates LDA and HDA at low-T, and continues into the liquid phase at high-T sep- 
arating the two liquid phases, LDL and HDL. Furthermore, it has been hypothesized that 
this first order transition line ends in a second (metastable) critical point [17,18,20,24,54]. 
The key points in the phase diagram proposed in [4,17,18] are that LDA is the glass obtained 
upon cooling of LDL while HDA is the glass obtained upon cooling of HDL. Recent results 
and simulations call for a refinement of this picture. In [30], we show that VHDA, and 
not HDA, is the glass obtained upon cooling of HDL. In this modified phase diagram, the 
first order transition line between LDL and HDL is continued in the glassy phase separating 
LDA and VHDA. We note that the properties of the glass obtained upon cooling a liquid 
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can change on the cooling rate [46]. Therefore, VHDA and LDA represent each a family of 
glasses, all of them obtained upon cooling (with different cooling rates) of HDL and LDL, 
respectively. Other unrelaxed glasses, e.g. HDA (and probably many other glasses [34] 
which might depend, for instance, on aging) cannot be shown in a single phase diagram. 
The phase diagram in [4,17,18] with a first order transition line identifies two distinct liquid 
phases and the two families of glasses associated to them. 
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TABLES 



Label 


P [g/cm 3 ] 


P [GPa] 


A 


0.90 


-0.55 


B 


0.94 


-0.17 


C 


0.96 


0.01 


D 


1.005 


0.41 


E 


1.065 


0.68 


F 


1.17 


0.84 


G 


1.255 


1.10 


H 


1.311 


1.38 


I 


1.375 


1.90 


J 


1.452 


2.66 


K 


1.532 


3.96 



TABLE I. Location in the P-p plane of the starting configurations used during the isobaric 
heatings. State points A-D correspond to LDA configurations while state points H-K represent 
HDA configurations. State points E-G correspond to configurations in the LDA-HDA transforma- 
tion. 
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FIGURES 




0.85 0.95 1.05 1.15 1.25 1.35 1.45 

p [g/cm 3 ] 

FIG. 1. Transformation of LDA to HDA obtained by isothermal compression at T = 0, 77, 120, 
and 170 K. LDA was obtained by hyper quenching at constant p equilibrium liquid configurations 
at T = 220 K and p = 0.9 g/cm 3 . The transformation starts approximately at (p c = 1-03 g/cm 3 , 
P c = 0.7 GPa) at T = K and shift with T to (p c = 0.97 g/cm 3 , P c = 0.13 GPa) at T = 170 K. 
We also show data for equilibrium liquid at T = 210 K (from Ref. [17]), one of the lower T 
accessible for equilibrium liquid simulations. We note that glass isotherms cross at (p ~ 0.98 g/cm 3 , 
P ~ 0.2 GPa) and that as T increases, the glass isotherms approach the liquid isotherm. 



27 



3.9 
3.4 



„ 2.9 

CL 

£ 2.4 

CL 



T=77K * 
T=170K 
T=210K 
T=230K 



Equil. 
Liquid 



1.9 
1.4 




0.9 



1.30 1.35 



1.40 



1.45 

p [g/cm 3 ] 



1.50 1.55 



FIG. 2. Magnification of Fig. 1 corresponding to densities of HDA. To test the presence of 
VHDA we heat at constant P configurations corresponding to state points H, I, J and K in the 
figure (see also Table I). Upon heating from T = 77 K (squares) up to T ~ 170 K (diamonds), 
the four state points move to higher densities crossing the glass isotherm for T = 170 K and 
approaching the T = 210 K liquid isotherm (see arrows). 
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FIG. 3. Evolution of p upon isobaric heating for the four state points H, I, J, K corresponding 
to HDA, indicated in Fig. 2. In all cases, p increases with T before the liquid phase is reached. The 
denser glass to which HDA transforms is identified as the VHDA obtained in experiments. Arrows 
indicate the increase in p when heating the glass from T = 77 K up to T 170 K. The increase in 
p is smaller as P increases. 
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FIG. 4. Effect of aging configurations of HDA corresponding to state point H in Fig. 2. With 
time, the p increases approaching the p corresponding to VHDA at P = 1.38 GPa. This finding is 
consistent with the view of VHDA as the result of relaxation of HDA upon isobaric heating. 
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FIG. 5. Magnification of Fig. 1 for intermediate densities in the LDA-HDA transition. We heat 
isobarically glass configurations corresponding to amorphous ices obtained during the transforma- 
tion of LDA to HDA (state points D, E, F), and one configuration corresponding to LDA (state 
point G; see also Table I). All state points shift to higher p upon heating (see arrows). Squares 
indicate the p of the starting glasses at T = 77 K while diamonds indicate the final p reached at 
T ~ 170 K before the glass transition. We also show the P — p values corresponding to two liquid 
isotherms (from Ref. [17]). 
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FIG. 6. Evolution of p upon isobaric heating for the four state points D, E, F, and G indicated 
in Fig. 5. Arrows indicate the increase in p when heating the glass from T = 77 K up to T ss 170 K. 
The smallest increase in p occurs for the heating of configurations indicated by state point G and 
it is larger for heatings at lower Ps. 
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FIG. 7. Magnification of Fig. 1 corresponding to densities of LDA. We heat at constant P 
configurations corresponding to state points A, B, and C in the figure (see also Table I). We show 
all the glass isotherms simulated corresponding to T = 0, 77, 100, 120, 140, and 170K (see dashed 
arrow), and the liquid isotherms for T = 210, 230, and 270 K (from Ref. [17]). Note that glass 
isotherms cross each other at (p ~ 0.98 g/cm 3 , P ~ 0.2 GPa). After heating and before the liquid 
phase is reached, state points A, B, and C shift to lower densities. In the three cases, the state 
points finally approach the liquid isotherm at high T. 
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FIG. 8. Evolution of p upon isobaric heating for the three state points A, B, and C indicated in 
Fig. 7. (a) Upon heating at P = 0.01 GPa, p decreases from 0.96 g/cm 3 at 77 K to p 0.952 g/cm 3 
at T = 150 K (see arrow) and state point C moves through the glass isotherm in Fig. 7. Then p 
increases and state point C shifts back to higher densities approaching the T = 230 K equilibrium 
liquid isotherm, (b) Upon heating the glass up to T « 170 K at P = —0.17 GPa, p decreases down 
to k, 0.93 g/cm 3 (see arrow). Upon farther heating up to 230 K, p does not change much probably 
because the p of the equilibrium liquid at P = —0.17 GPa and 230 K is also ss 0.93 g/cm 3 . (c) 
In this case, p decreases monotonically and state point A in Fig. 7 shifts to lower densities. The p 
reaches very low values and the glass transforms into gas. Arrow indicates the increase in p when 
heating the glass from T = 77 K up to T w 160 K. 
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FIG. 9. (a) P during decompression of the system starting at state points H, I, H" and I". Inset: 
State points H' and I' indicate VHDA configurations at T ~ 165 K obtained after heating the HDA 
configurations represented by state points H and I at T = 77 K. State points H" and I" correspond 
to the final VHDA states after cooling VHDA back to 77 K. Decompressions at T = 77 K of (b) the 
four HDA states indicated by state points (H, I, J, K) and (c) the four VHDA states indicated by 
(H", I", J", K"). At P « GPa, the recovered configurations of HDA and VHDA do not collapse 
into two single states but they cover two ranges of densities. HDA densities fall in the interval 
1.15 — 1.24 g/cm 3 while VHDA densities fall in the interval 1.22 — 1.28 g/cm 3 . All recovered HDA 
and VHDA forms collapse to a single state at p = 1.05 g/cm 3 and P ~ —0.4 GPa. 
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FIG. 10. Oxygen-oxygen (OO), oxygen-hydrogen (OH), and hydrogen-hydrogen (HH) radial 
distribution functions at = 77 K and P = 0.02 GPa for (a) LDA and HDA, and (b) HDA and 
VHDA. The area under the RDFs indicate that LDA, HDA and VHDA are characterized by 
a tetrahedral hydrogen-bonded network and that, in comparison with LDA, HDA has an extra 
interstitial neighbor between the first and second shell while VHDA has two extra molecules. 
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FIG. 11. Evolution of the radial distribution functions upon isobaric heating of 
HDA(1.37 g/cm 3 ) at P = 1.9 GPa from T = 77 K up to T « 175 K. 
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FIG. 12. (a) Compression of LDA to obtain HDA, and subsequent decompression of HDA from 



1.4 g/cm 3 down to 0.85 g/cm 3 . Triangles indicated the starting configurations of recovered HDA 



at (1.2 g/cm 3 , 0.01 GPa) and (1.08 g/cm 3 , 0.01 GPa) that are heated at constant P. The p of 



the system upon heating is shown for (b) P = 0.01 GPa and (c) P = —0.3 GPa. In both cases, 



the system evolves toward the T = 210 K liquid isotherm. Arrows in (a) indicate the p change 



upon heating up to T = 170 K, before the liquid phase is reached. Only the isobaric heating at 



P = —0.3 GPa allows the system to reach the p corresponding to LDA (~ 0.93 g/cm 3 for this P) 



supporting the possibility of a HDA— s-LDA transformation. 
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FIG. 13. (a) LDA (p = 0.9 g/cm 3 ) is compressed at T = 77 K to obtain HDA (p » 1.31 g/cm 3 ). 
Then, HDA is heated isobarically (horizontal dotted arrow) to obtain VHDA (filled triangle down) 
and expanded at T = 77 K (dashed line). VHDA recovered at p = 1.27 g/cm 3 (upper filled 
triangle left) is heated at constant P = 0.11 GPa to test whether it produces LDA (see panel 
(b)). At this P, the system gets trapped in the liquid phase and no transformation to LDA is 
observed. VHDA recovered at p = 1.11 g/cm 3 (lower filled triangle left) is heated at constant 
P = —0.3 GPa. The evolution of p upon heating (see panel (c)) is consistent with the possibility of 
a VHDA^LDA transformation. VHDA recovered at p = 1.26 g/cm 3 (filled triangle up) is heated 
at constant p = 1.26 GPa to test whether it produces HDA (see panel (d)). The location of the 
T = 210 K liquid isotherm avoids a VHDA^HDA transformation (at the present heating rate 
q h = 3x 10 10 K/s). 
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